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Abstract 

Background: In longitudinal studies where subjects experience recurrent incidents over a period of time, such as 
respiratory infections, fever or diarrhea, statistical methods are required to take into account the within-subject 
correlation. 

Methods: For repeated events data with censored failure, the independent increment (AG), marginal (WLW) and 
conditional (PWP) models are three multiple failure models that generalize Cox's proportional hazard model. In this 
paper, we revise the efficiency, accuracy and robustness of all three models under simulated scenarios with varying 
degrees of within-subject correlation, censoring levels, maximum number of possible recurrences and sample size. 
We also study the methods performance on a real dataset from a cohort study with bronchial obstruction. 

Results: We find substantial differences between methods and there is not an optimal method. AG and PWP seem to 
be preferable to WLW for low correlation levels but the situation reverts for high correlations. 

Conclusions: All methods are stable in front of censoring, worsen with increasing recurrence levels and share a bias 
problem which, among other consequences, makes asymptotic normal confidence intervals not fully reliable, 
although they are well developed theoretically. 



Background 

Coxs proportional hazards model [1] is the most com- 
monly used model for clinical trial data and provides 
reliable estimates of survival times, as well as the rela- 
tive risk associated with time-to-event occurrence. As a 
semi-parametric model, it does not have any constraints 
on distributional assumptions, which makes it an attrac- 
tive alternative to parametric models. However, it is not 
suitable for recurrent event data because survival times in 
the standard model terminate at the time of the event. In 
addition, trying to use multiple failure times is inappropri- 
ate since events within individuals may be correlated. In 
this case, the assumption of independence in the standard 
Cox regression model is violated, introducing statistical 
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complications. To avoid the errors resulting from ana- 
lyzing correlated repeated events, sometimes repeated 
events are disregarded and time to first event is only 
used, which is clearly a suboptimal approach. Multiple- 
event analytical techniques for survival data have been 
developed during the last decades [2] but, as a result of 
their complex structure and computational requirements, 
they have not been commonly applied. Advancement in 
statistical software in recent years has made these meth- 
ods more accessible to researchers. As a consequence, 
their use has gradually increased in areas such as clin- 
ical research. Text books like Therneau and Grambsch, 
[2], describe these models in detail with examples, while 
research articles, such as that of Kelly and Lim [3], pro- 
vide a very comprehensive explanation of their application 
to multiple-failure data, with extensive discussion about 
assumptions. There are three commonly used statistical 



O© 201 3 Villegas et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative 
BlolVted Cental Commons Attribution License (http://creativecommons.Org/licenses/by/2.0), which permits unrestricted use, distribution, and 
reproduction in any medium, provided the original work is properly cited. 



Villegas et al. BMC Medical Research Methodology 201 3, 1 3:95 
http://www.biomedcentral.com/1471-2288/13/95 



Page 2 of 10 



methods for the analysis of recurrent events: Anderson- 
Gill model [4], Prentice-William-Peterson model [5] and 
Wei-Lin and Weissfeld model [6]. However, the benefits 
and drawbacks of each method are still unclear to the 
medical researchers, specially under varying strengths of 
correlation between survival times. The aim of this work 
is to study, through simulations and a real data set, the dif- 
ferences between the three models under presence of low 
to moderate correlation in survival times. With respect to 
similar studies like Kelly and Lim [3], we very approxi- 
mately control the within individual correlation between 
consecutive recurrent events and we allow decreasing cor- 
relations as long as distance between recurrent events 
increases. In addition our simulation design takes care of 
different censoring levels and different numbers of recur- 
rences, which are not fixed at 4 events as in [3]. We briefly 
describe the models, discuss in more detail the simulation 
procedure, and discuss the results of fitting the simulated 
data for the different models. In the Example section, we 
provide an illustration of the models performance using 
real data from a study in infants with respiratory illness 
recurrences. 

Methods 

Basic notations 

It is common to model survival times through the haz- 
ard function. The Cox proportional hazard model for right 
censored survival data specifies the hazard function for 
the i th individual by 

hi(t) = ho(t)exp(fi'Zi), (1) 

where t represents time, Z; is the vector of covariates of 
the i th individual, ft is a vector of regression coefficients 
and ho(t) is the so-called baseline hazard function that 
corresponds to the hazard function for an individual with 
covariables set to zero. As model (1) is formulated through 
the hazard function, the simulation of appropriate survival 
times for this model is not straightforward, the effect of 
the covariates must be translated from the hazards to the 
survival time. We need to simulate the survival times since 
software packages for Cox models require the individ- 
ual survival time data, not the hazard function. Suppose 
that there are n individuals and that each individual can 
experience K failures. Times of failure are subject to right 
censoring and we consider that censoring is non informa- 
tive, which means that knowledge of a censoring time for 
an individual or subject provides no further information 
about the individuals probability of survival at a future 
time. Let Ty be the time when the j th failure occurs for the 
i th subject, measured from the subjects study enrollment, 
and Q be the corresponding censoring time for individual 
i. Since the j th failure time is affected by right censoring 
we define Xy as the minimum of (7^, Q), i.e., Xy equals 
Tu if the event was observed, and Q if it is censored. Let 



8 t j = I(Tij < Q) be the indicator function taking the value 
1 when Tij < Ci and 0 otherwise. For simplicity, in our 
simulations all censoring times are taken as constant, i.e. 
Q = C for all /. 

Overview of the multiple failure time models 

For multiple events the common approaches that gener- 
alize the Cox's framework are the independent increment 
(AG), marginal (WLW) and conditional (PWP) models. 
These three models differ essentially in the risk sets. 

Andersen and Gill (AG) independent increment model 

In the AG model [4] it is assumed that recurrent events 
are unaffected by earlier events that occurred to the same 
subject so baseline hazards for all events are common. The 
assumption of mutual independence of the events within 
a subject is equivalent to the assumption of independent 
increments in the counting process inside each individ- 
ual Each recurrent event for the i th subject is assumed 
to follow a proportional hazard model where the hazard 
function is 

hi(t) = ho(t)exp(ft'Zi(t)). 

Under this model, the risk of a recurrent event for a sub- 
ject follows the usual proportional hazards assumption 
but the rank of recurrence is not taken into account. By 
defining the risk set indicator Yy(t) = liX^-i < t < Xij), 
the risk set at time £, represented by Yij(t), includes 
all subjects under observation regardless of the number of 
occurrences experienced by each subject. As j increases, 
fewer individuals are likely to experience the event, but 
this fact does not affect the risk set due to the fact that it 
does not depend on so the coefficients estimates are sta- 
ble. The AG model provides more efficient inference for 
a covariate effect than the plain Cox model for the time 
to the first event, but requires much stronger assump- 
tions than the other models, such as independence among 
recurrent event times or common baseline hazards. Also, 
a robust "sandwich" method can be used in the estimation 
of standard errors [7]. 

Prentice, Williams and Peterson (PWP) conditional model 

Another model for analyzing recurrent events is the PWP 
model [5] with time scale measured from the beginning of 
the study to a specified failure. It assumes that a subject 
is not at risk for the j th event until he/she has experienced 
event / — 1. This model requires the same assumptions as 
the Cox model, but, unlike the AG model, it allows the 
baseline hazard to vary from recurrence to recurrence, the 
hazard function for the j th event for the i th subject is: 

hij(t) = hoj(t)exp(PjZi(t)). 

The risk set indicator is the same as in the AG model, 
Yij(t) = I(Xij-i < t < Xij), but the risk set at time t 
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is different for each j: ^ According to the defini- 

tion of the risk set, the number of subjects is dramatically 
decreased as j increases. Stable coefficient estimates can- 
not be obtained for higher ranks of The hazard function 
at time t for the j th recurrence is conditional to the entire 
previous failures. This model allows different baseline 
hazards, therefore, estimations for the current event may 
be affected by earlier events. 

Wei, Lin, and Weissfeld (WLW) marginal model 

The WLW model [6] uses, for each the distribution of 
time from the study enrollment to the f h recurrent event 
according to a Cox model, leaving the dependence struc- 
ture between the observations within the same subject 
completely unspecified. For the WLW model the risk set 
indicator is Yiy(£) = I(Xij > t), and the risk set at time t 
for the j th recurrence, represented by Yyit), includes 
anyone who has not experienced the j th recurrence at time 
t. In other words, unlike the PWP model, a subject in 
the risk set at time t does not necessarily have to experi- 
ence the (j — l) th recurrence. The hazard function of the 
marginal model for the j th event for the i th subject is 

hijit) = h 0j (t)exp(PjZi(t)). 

The WLW approach obtains the estimate of the covari- 
ate effects from the partial likelihood function under the 
working assumption of independence and then adjusts the 
variance estimate empirically. Although observations on 
each subjects are correlated, the ft estimation has been 
shown to be consistent in the case of two events [8]. As 
the inverse of the information matrix provides a naive 
variance-covariance estimation, not correctly estimating 
the true variance, a robust "sandwich" variance-covariance 
estimation method is preferred. 

In all these models we can add a frailty, that is, a random 
individual effect to explain the dependence of individ- 
ual recurrent events. In our simulation study we try to 
precisely control the degree of correlation between adja- 
cent recurrent event times. Frailty makes difficult such a 
control, so we have not considered it. 

Simulation 

We carried out a series of simulations to examine the 
accuracy of the above multiple failure time models in 
terms of four factors: different sample sizes, censoring 
levels, numbers of recurrence events and correlation lev- 
els. The sample sizes under consideration were n — 
(50, 100,200,400). Survival times were censored at fixed 
time determined to yield censoring percentages p of 0%, 
15%, 30% and 50%. In other words, previously fixed cen- 
soring times C p = (Co = +oo, Q5, C30, C50) were estab- 
lished in order to ensure the preceding proportions of 
subjects experiencing a censure. The maximum num- 
ber of recurrent events under consideration were K = 



(3, 6, 9, 12). For subject /, the recurrence times tn, . 

were generated as correlated Weibull deviates as is 
explained below. Censoring occurred in subject i if tn + 
fe + • • • + Uk > C p , In this case he/she experienced / 
events, j < K, if Ty = + • • • + < C p and tn + . . . + 
> C p . Note that not always censoring occurred only 
in the K th event. The correlation levels between adjacent 
recurrence times were set at p = (0, 0.10, 0.45, 0.80). The 
320 simulated scenarios were defined by crossing all lev- 
els of these four factors. Each simulation series consisted 
of the generation of m = 10000 random samples, each 
one representing a survival dataset where half of individ- 
uals were considered as "control" and half of individuals 
were considered as "treatment". This covariate was rep- 
resented by a binary variable with values 0 and 1. The 
treatment effect was set at ft = —1. For each simulated 
dataset we adjusted the WLW, AG and PWP models. We 
used the following measures of the accuracy of the regres- 
sion parameter estimates and the gain of robustness from 
the robust approaches: the true coverage of the 95% con- 
fidence interval for the parameter ft, the relative bias of 
the /3 estimators under the three models and the bias of 
the "robust" and "naive" variance estimators for all three 
models. The relative sampling bias is defined as the aver- 
age bias from the m random samples. That is, if fa is the 
estimate on i th simulated dataset, then 

1 m (fa fa) 

relative sampling bias = — 7 1 



i=l 



To simulate the survival times that meet the proportional 
hazards assumption we used the algorithm proposed by 
Feiveson [9]. Let U be a random variable with uniform dis- 
tribution. We can generate a new random variable Y with 
distribution function F by using the inverse (or pseudo 
inverse) function F~ l [10] as follows 

Y = F~ 1 (U). 

In our case, we are interested in simulating times with 
Weibull distribution with parameters c and k, that is, with 
distribution function F(y) = 1 — exp (— cy K ), for c > 0, 
k > 0. Note that, usually, the Weibull parameters are X 
and k with c = X K . Therefore, given U ~ Uniform [0, 1], 



(2) 



has the aforementioned Weibull distribution. The hazard 
function corresponding to Y is given by h(y) = cicy K ~ l . 
Under the Cox model and for a binary covariate X, the 
parameter c would be equal to some base value cq when 
X = 0 and would become equal to c\ = cq exp(^) when 
X = 1. 
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Simulation of within-subject dependence 

Our goal is to simulate recurrent Weibull failure times 
with correlation p between successive recurrence events 
in the same individual and with decreasing correlation 
as the distance between events increases. To incorpo- 
rate this dependence into our simulated data while still 
maintaining the marginal Weibull distribution for each 
observation we will use a transformation of a convenient 
multivariate normal distribution. We describe this proce- 
dure in three steps: 

1. For any correlation p, choose an appropriate po fol- 
lowing formula (9) derived in the next section. For 
each subject, i, i = 1, . . . , n, generate a random vari- 
ate Zn from a N(0, 1) distribution. Given zn, generate 
a random variate Zq from a normal distribution with 
mean poZn and variance 1 — p%. Given za, gener- 
ate a random variate z^ from a N(poZi2, 1 — p%); 
and so on, until zik> This ensures that each con- 
secutive pair (zij,Zij+i) is a realization of (Z^Z^+i) 
following a bivariate normal distribution with stan- 
dard marginals, N(0, 1) and correlation po. Also the 
correlation between Z,y and Zj,y+ T decreases as r 
increases. 

2. Transform Zy into uniform random variables using 
the standard normal cumulative distribution func- 
tion O: Uij = O (Zij). Variables Uij have a similar 
correlation structure. 

3. Finally, transform the variables Uij into dependent 
Weibulls using Equation (2). These Weibulls have 
again a similar correlation structure, with correlation 
p between adjacent recurrence times. 

Computation of the correlation p 0 

In order to simplify our argumentation we use the loga- 
rithm of Weibull times. Using delta method it is easy to 
see that the correlation between times is approximately 
equal to the correlation between time logarithms. In terms 
of the original normally distributed Z,y, the simulated log 
failure times may be expressed asg// = G (Z^), where 



G(z) 



r-iog<i>(z)"|if 



(3) 



Using second-order Taylor expansions of and gij+i 
about z = 0, it can be shown (see Appendix) that 



Cov ( gij ,gi,j +1 ) « {G'(0)} 2 po + ^ ^ Pi (4) 



and 



Var(g lj ) = Var(gi, j+l ) « {G'(0)} 2 + 



{G"(0)} 2 



(5) 



Combining (4) and (5), we find that Corr (gij>gi,j+i) is a 
weighted average of po and p^\ namely, 



Corr (gij,gi,j+i) « wp 0 + -(1 - w)p% 



where 



{G'(0)Y 



{G'(0)Y + {G"(0)} z /2 



(6) 



(7) 



It can be shown (see Appendix) that w does not depend on 
the parameters k or c; specifically, 



Ti (log2) 2 



^(log2) 2 + (l-log2) 2 



: 0.94128 



(8) 



Thus, solving (6) for the value of po such that 
Corr (gij,gi,j+i) = p, the solution is 



Po 



-w + yV 2 + 2p (1 - w) 
(l-w) ' 



(9) 



Software simulation 

All the simulations were performed using the library sur- 
vival and the functions qnorm and runif from the sta- 
tistical software R2.15.2, 64 bit version, [11]. The default 
generator is based on a Mersenne-Twister algorithm. 

Results 

Table 1 shows results for a simple size of 200 subjects and 
15% censoring. For simplicity, we provide only graphical 
results for two degrees of correlation between adjacent 
recurrent event times (0.1 and 0.8) and two maximum 
numbers of possible recurrent events (3 and 9). The com- 
plete simulation results corresponding to all 320 simula- 
tion scenarios are accessible at the url http://www.ub.edu/ 
stat/recerca/materials/SupplEmpirical.htm. 

If only the first event is considered, the relative bias, 
efficiency and coverage of the standard Cox proportional 
hazard estimator are adequate, as expected, in all the 
investigated scenarios. We do not report these results 
because our interest is focused in recurrent data. 

Censoring has a very small influence on the perfor- 
mance of the treatment effect estimators, /3, for all models 
under consideration. The graphics are nearly identical for 
all levels of censoring so only 15% level is reported. 

Figure 1 illustrates the relative bias, Figure 2 the bias 
of the variance estimators and Figure 3 the true coverage 
of the 95% confidence interval based on the robust vari- 
ance estimator. All these figures are organized in rows and 
columns. Each column corresponds to a correlation level 
(p = 0.1 and p = 0.8) and each row to a maximum 
number of possible recurrences (K = 3 and K = 9). 

With respect to relative bias, the AG estimator of /3 is 
stable for all different sample sizes. The bias grows with 
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Table 1 Relative bias, naive variance, robust variance and coverage robust for 1 5% censoring and sample size of 200 

p Events Model Relative bias Naive var Robust var Coverage 

0.1 3 AG -0.193 0.008 0.006 0.308 

PWP -0.030 0.010 0.011 0.930 

WLW 0.320 0.009 0.021 0.415 

6 AG -0.254 0.004 0.003 0.003 

PWP -0.042 0.006 0.007 0.897 

WLW 0.666 0.005 0.021 0.004 

9 AG -0.278 0.003 0.002 0.000 

PWP -0.048 0.005 0.006 0.871 

WLW 0.929 0.004 0.022 0.000 

12 AG -0.291 0.002 0.001 0.000 

PWP -0.050 0.004 0.005 0.861 

WLW 1.151 0.003 0.022 0.000 

0.45 3 AG -0.296 0.008 0.007 0.069 

PWP -0.201 0.009 0.011 0.491 

WLW 0.166 0.009 0.022 0.806 

6 AG -0.343 0.004 0.004 0.001 

PWP -0.291 0.005 0.006 0.055 

WLW 0.358 0.005 0.022 0.318 

9 AG -0.350 0.003 0.003 0.000 

PWP -0.324 0.004 0.004 0.008 

WLW 0.514 0.003 0.021 0.058 

12 AG -0.349 0.002 0.002 0.000 

PWP -0.338 0.003 0.003 0.002 

WLW 0.657 0.003 0.022 0.005 

0.8 3 AG -0.381 0.008 0.008 0.016 

PWP -0.339 0.009 0.011 0.124 

WLW 0.051 0.009 0.024 0.936 

6 AG -0.448 0.004 0.005 0.000 

PWP -0.509 0.005 0.005 0.000 

WLW 0.110 0.004 0.023 0.892 

9 AG -0.462 0.003 0.004 0.000 

PWP -0.574 0.003 0.003 0.000 

WLW 0.157 0.003 0.023 0.829 

12 AG -0.462 0.002 0.004 0.000 

PWP -0.605 0.002 0.003 0.000 

WLW 0.205 0.002 0.022 0.727 



the correlation and the number of possible recurrences. 
The treatment effect is always underestimated. In all cir- 
cumstances the naive variance estimator is worse than the 
robust estimator, but both perform very similarly and ade- 
quately. The true coverage of the confidence intervals is 
in general low and decreases with correlation, number of 
events and, surprisingly, with sample size. 



The estimator of f$ based on model PWP is nearly 
unbiased for low correlations and stable with respect 
to the number of recurrences. When the correlation is 
high, it underestimates the treatment effect /3. This bias 
does not improve with growing sample sizes. Again, the 
robust variance is better than the naive variance, but both 
are very similar and perform adequately, improving with 
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Relative bias for the AG, PWP and WLW model 

with p=0.1, censoring=15% and 3 events 



Relative bias for the AG, PWP and WLW model 

with p=0.8, censoring=15% and 3 events 
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Relative bias for the AG, PWP and WLW model 

with p=0.1, censoring=15% and 9 events 



Relative bias for the AG, PWP and WLW model 

with p=0.8, censoring=15% and 9 events 
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Figure 1 Relative Bias for each model with correlation of 0.1 and 0.8, 1 5% censoring and 3 and 9 events. 



growing sample sizes and correlations. For low correla- 
tion levels its coverage is the best, but it quickly falls with 
increasing correlations and becomes sensitive to growing 
sample sizes, in the same surprising direction as the AG 
method. 

Under the WLW model is overestimated. This bias 
increases with growing recurrence but improves with 
growing correlation. The naive variance estimator under 
WLW is very biased, its use should be avoided. The robust 
variance estimator is clearly better but slightly worse than 
the other robust estimators for sample sizes below 200. 
Under low correlation, the confidence interval performs 
very inadequately, with very low coverages which decrease 
with growing sample sizes and with growing number 
of recurrences. On the other hand, its coverage greatly 
improves with growing correlations, though a tendency 
to make worse with growing sample sizes and recurrence 
levels still persists. 

Example 

The models under consideration are illustrated through 
the analysis of a cohort of infants with bronchial 



obstruction [12]. For this study, 4 months-old children 
were followed until they reached 12 months of age and 
recruitment took place during the course of 18 months, 
from April 1995 to October 1996. The chief goal of the 
investigation was to estimate the effects of fine particu- 
late matter and wheezing illness in the first year of life. 
The recurrent events of interest were physician office 
visits attributable to respiratory illnesses. From the pedi- 
atric visit database, a total of 504 infants were eligible 
for the present study in the pediatric visit database. The 
events of interest were the times to physician revisit fol- 
lowing the initial visit due to examination. Each subject 
experienced a particular number of visits to physician at 
various times, which represent the whole observable his- 
tory of his/her recurrences. For each individual the last 
time is censored due to the end of the study period. Of 
these 504 subjects, 475 had at least one revisit during the 
follow-up time. Table 2 summarizes the number of events 
experienced by the 475 subjects during the follow-up 
time. A total of 475 occurrences on revisit were observed, 
where some infants experienced a rather large number 
(up to 10 revisits). 
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Estimated bias for naive and robust variances 
in the AG, PWP and WLW model 



with p=0.1, censoring=15% and 3 events 



Estimated bias for naive and robust variances 
in the AG, PWP and WLW model 



with p=0.8, censoring=15% and 3 events 




Estimated bias for naive and robust variances 
in the AG, PWP and WLW model 

with p=0.1, censoring=15% and 9 events 



Estimated bias for naive and robust variances 
in the AG, PWP and WLW model 

with p=0.8, censoring=15% and 9 events 
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Figure 2 Estimated Bias for naive and robust variances for each model with correlation of 0.1 and 0.8 7 1 5% censoring and 3 and 9 events. 



The main purpose of this example is to illustrate 
the application of recurrent event techniques and pro- 
vide comparisons between techniques, rather than giv- 
ing estimates for recurrence of respiratory illnesses. Two 
characteristics, parents reported smoking (l=yes, 0=no) 
and family history of asthma (l=yes, 0=no) were included 
as two univariate analysis in each of the models studied. 
We chose these two covariates because they are poten- 
tially related with the outcome as indicated in the medical 
literature and also because they are binary variables as we 
used in our simulation study. Table 3 contains the results 
of three different models: AG, PWP and WLW, fitted to 
the respiratory data. The /3 values estimated by the AG 
and PWP models are lower than those estimated by the 
WLW model for both covariates, smoking and asthma. 
But all of them show a negative effect in the outcome. 
When we analyze the standard errors, the robust standard 
errors are larger in all models than those estimated by the 
naive method. And again, the WLW model has the biggest 
standard error in comparison to the other models. These 
results are coherent with the previous simulation results. 



Note that the WLW is the only model where asthma is a 
statistically significant. 

Discussion 

The simulation results show that each model has dif- 
ferent degrees of robustness against variations in the 
number of recurrent events and correlation. In terms of 
relative bias, the WLW model shows the worst perfor- 
mance for low correlations but greatly improves for high 
correlation levels, independently of censoring and sam- 
ple size. The situation is reversed for the AG and PWP 
models, which perform better than WLW for low correla- 
tion levels. However, these models clearly underestimate 
the parameter. As some authors have suggested [3], the 
WLW model overestimates regression coefficients due to 
a carry-over effect. This is corroborated by our simu- 
lation results. Additionally, as we can see from Table 3 
in the example, the WLW model presents the highest 
coefficients and robust standard errors. 

The most surprising results refer to the unexpected 
behavior of all methods with growing sample sizes. First, 
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Coverage robust 95% CI for the AG, PWP and WLW model Coverage robust 95% CI for the AG, PWP and WLW model 

with p=0.1 , censoring=15% and 3 events with p=0.8, censoring=15% and 3 events 
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Coverage robust 95% CI for the AG, PWP and WLW model Coverage robust 95% CI for the AG, PWP and WLW model 

with p=0.1 , censoring=15% and 9 events with p=0.8, censoring=15% and 9 events 
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Figure 3 Coverage robust 95% confidence intervals for each model with correlation of 0.1 and 0.8, 15% censoring and 3 and 9 events. 



their bias remains nearly constant. This was so surpris- 
ing that we repeated the computations with Stata, with 
identical results to those obtained with the "coxph" R 
function. Provided that (as is expected) the true vari- 
ance of all estimators decreases and its robust estimation 
greatly improves with growing sample sizes (Figure 2), as 
a consequence the confidence intervals become more and 
more narrow around a biased value, which produces a 
disconcerting, at first sight, fall of coverage (Figure 3). 

The AG, PWP and WLW models have been previously 
compared by various authors (e.g., [3,7,13,14]) using real 
and simulated data, showing that these models often 
yield different results for the same data set. This is 
not unexpected, since distinct models address different 
research questions. For example, in [3] they found that 
the WLW model overestimates the treatment effect. Also 



they show that as the correlation increases the coverage of 
the 95% confidence interval based on the robust variance 
decreases, except for WLW. 

Conclusion 

In this paper we reviewed and compared the principal 
methods in the analysis of recurrent event data in 
terms of censoring and within-subject correlation. The 
example illustrates the differences between the methods 
and resulting parameter estimators. 

No method can be recommended as the best in all 
circumstances. The AG and PWP approaches are quite 
robust in front of low levels of within-subject correlation, 
but the WLW approach should be recommended under 
the suspicion of high correlation. All methods share a 
bias problem which makes the confidence intervals based 



Table 2 Number of episodes of illness 



Number of events 




0 


1 


2 


3 4 


5 


6-10 


Total 


N. of subjects (%) 


260 (54.7) 


101 (21.3) 


53(11.2) 


21 (4.4) 21 (4.4) 


10(2.1) 


9(1.9) 


475 (100) 
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Table 3 Coefficients and standard errors obtained in two univariate analysis for 2 covariates: parents reported smoking 
and family history of asthma 



Model 


Covariate 


Beta 


Naive S.E. 


Robust S.E. 


Robust 95% CI 


AG 


Smoking 


0.158 


0.059 


0.089 


-0.016,0.332 




Asthma 


0.276 


0.111 


0.155 


-0.028 , 0.579 


PWP 


Smoking 


0.147 


0.060 


0.088 


-0.025,0.319 




Asthma 


0.251 


0.113 


0.148 


-0.039 , 0.540 


WLW 


Smoking 


0.162 


0.058 


0.102 


-0.037 , 0.362 




Asthma 


0.407 


0.111 


0.189 


0.036 , 0.778 



on asymptotic normal theory not fully accurate. On the 
other hand, the robust estimation of the variance of the 
treatment effect estimator is very reliable, and should be 
recommended instead of the naive variance estimation, 
especially for the WLW method. 

Appendix 

Remember that p{Zij,Z ik ) = p 0 > 0 and p(gij,gi,j+i) = 
p > 0 with gij = G(Zij) where G is defined in (3). To 
find the correlation between gij and gi,j+\ we will use a 
second order Taylor expansion (delta method, [15]) for the 
function G about z = 0 as follows: 



G"(0) 



gy = G(Zy) = G(0) + G f mZij - 0) + ^^(Z, y -0) 2 + . 
G(0) + G f (0)Zij + ^y^Z? 



G"(0) , 



since Z,y ~ N(0, 1), we have E(Zij) = 0 and VariZij) = 



E(Zp = 1, then, 



Eigij) « G(0) + 



G"(0) 



gij ~ Eigij) « G\0)Zij + ^^(Z? - 1) 
g Uj+l - E(g Uj+l ) » G'(0)Z,* + ^^(Zf k - 1) 



and 



igij-E{gij))(g Uj+ i-E{g iij+l )) « (G\0)yZijZ ik 

G f (0)G ff (0) 0 G'(0)G"(0) o 



+ 



2 

(G"(0)) 2 9 9 

^^(z 2 -i)(4-i). 



By applying now the £(•) operator taking into account that 



so £[ ZyZ tt ] = po and E[ Z tj {Z\ - 1)] = E[ Z ik Z]j] = 0, we 
have: 

9 (G"(0)) 2 
Cov^, &/+1 ) « (G'(0)) 2 £(ZyZ iVt ) + ' 

x£[(Z?-l)(4-l)] 

. 2 (G"(0)) 2 2 

= (G'(0)) 2 p 0 + ~ Pi 
For the Varfgy) we proceed in the same way as before, 

Varigij) « Vfcr (g(0) + G'(0)Z, 7 + ^y^Z 2 ) 

= (G\0)) 2 Var(Zij) + i(G>)) 2 Vfcr(Z?) 

= (G / (0)) 2 + ^(G // (0)) 2 . 

because Cov(Z iJf Zf k ) = 0, and Var(Z]j) = 2. Finally the 
formula for the approximated correlation between gij and 

gi,j+i is 



Corr(gij,g if j+i) = 



Cov(gij>gi,j+i) 
jVar(gij)Var(g i}j+ i) 

(G / (Q)) 2 po + |(G // (Q)) 2 p 0 2 
(G /(0))2 + i(G^(0)) 2 



2 V 
::Vrv^2 



(G'(0)^ 



(G'(0)) 2 + *(G"(0)) 2 



Po 



1 

+ - 



2 V 



(G"(0)r 2 



4(G'(0)) 2 + ±(G"(0)) 2 
= w • p 0 + — - — P 0 



where, w = (G/(q) |^T ( g// (0)) 2 ' In order to compute w, we 



use that: 



* (z) = f-00 -j^f hdx $(0) = 2 

= ^e-# ^ O'(0) = ^ 

<&"(*) = — |=e _ T O"(0) = 0. 

V 27T 
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Then, 



G'{z) 
G'(0) 



1 &(z) 
k <J>(z)log<J>(z) 

V2 



£(log2)V^ 

1 (<$>'(z)) 2 , 1 <$>"(z)<$>(z)-{<S>'{z)) 2 



w /< (log O (z)) 2 O (z) 2 £ O (z)2 log $ ( Z ) 



G"(0) 



And finally, 



2(log2-l) 

/T7T(log2) 2 ' 



(G'(0)) 2 



7T(l0g2) 2 



(G'(0)) 2 + ±(G"(0)) 2 7r(log2) 2 + (log2 - l) 2 ' 
Note that w does not depend neither on k nor on c. 
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